survivalパッケージのsurvfit関数は、生存時間データのカプランマイヤー推定で用いる関数である。 左側切断があるデータで、想定外の出力で困ったので、その対処法のメモ。

簡単な擬似データ

時点beginで観察開始（左側切断）、時点endで観察終了、観察終了時の状態がstatusである。

library(survival)

dat <- data.frame(
  id=0:5,               # ID
  begin=c(0,1,2,3,4,5), # 観察開始時点
  end=c(3,4,5,6,7,8),   # 観察終了時点
  status=c(1,0,1,1,1,1) # 観察終了時の状態（1：イベント発生、0：打ち切り）
)
dat
##   id begin end status
## 1  0     0   3      1
## 2  1     1   4      0
## 3  2     2   5      1
## 4  3     3   6      1
## 5  4     4   7      1
## 6  5     5   8      1

模式図

擬似データにおける追跡調査の様子の模式図。

時点0から観察しているとした場合（左側切断なし）

summary(km1)はカプランマイヤー推定による表である。

km1 <- survfit(Surv(end,status)~1,data=dat)
summary(km1)
## Call: survfit(formula = Surv(end, status) ~ 1, data = dat)
## 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     3      6       1    0.833   0.152       0.5827            1
##     5      4       1    0.625   0.213       0.3200            1
##     6      3       1    0.417   0.222       0.1468            1
##     7      2       1    0.208   0.184       0.0368            1
##     8      1       1    0.000     NaN           NA           NA

km1$timeでイベント発生時点、km1$survでその時点での生存率が抽出できると思っていたけど、違ったみたい。 結果から見れば、打ち切り時点である時点4が含まれているので、km1$timeは単にユニークな時間を出力するみたい。

km1$time
## [1] 3 4 5 6 7 8
km1$surv
## [1] 0.8333333 0.8333333 0.6250000 0.4166667 0.2083333 0.0000000

beginの時点から観察しているとした場合（左側切断あり）

summary(km2)はカプランマイヤー推定による表であり、ちゃんと左側切断の情報も考慮しているのが分かる。

km2 <- survfit(Surv(begin,end,status)~1,data=dat)
summary(km2)
## Call: survfit(formula = Surv(begin, end, status) ~ 1, data = dat)
## 
##  time n.risk n.event entered censored survival std.err lower 95% CI
##     3      3       1       1        0    0.667   0.272       0.2995
##     5      3       1       1        0    0.444   0.257       0.1433
##     6      3       1       0        0    0.296   0.210       0.0741
##     7      2       1       0        0    0.148   0.148       0.0209
##     8      1       1       0        0    0.000     NaN           NA
##  upper 95% CI
##             1
##             1
##             1
##             1
##            NA

左側切断のある生存時間データの場合、km2$timeで出力されるのは、左側切断時点と観察終了時点のユニークな時点を出力するみたい。 これに気付かずに詰んでた。

km2$time
## [1] 0 1 2 3 4 5 6 7 8
km2$surv
## [1] 1.0000000 1.0000000 1.0000000 0.6666667 0.6666667 0.4444444 0.2962963
## [8] 0.1481481 0.0000000

イベント発生時点のユニークな時点とその生存率が欲しい場合は、下記のようにすれば、とりあえず得られる。

s2 <- summary(km2)
s2$time
## [1] 3 5 6 7 8
s2$surv
## [1] 0.6666667 0.4444444 0.2962963 0.1481481 0.0000000